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1 Introduction 


We consider numerical methods for the time integration of a family of Initial 
Value Problems in ODEs 

y'hit) = fh{t,yh{t)), yh{0)=u*Qh, 0 <t<t*, yh,fheW^^’^\ 

( 1 . 1 ) 

coming from the spatial semi-discretization of an Z—dimensional Advection 
Diffusion Reaction problem in time dependent Partial Differential Equations 
(PDEs), with prescribed Boundary Conditions and an Initial Condition. Here 
h denotes a small positive parameter associated with the spatial resolution 
and usually I = 2,3,... . 

The typical PDE problem with Dirichlet boundary conditions is given by (D is 
a bounded open connected region in MJ, dfl its boundary and V is the gradient 
operator) 

Ut{x, t) = —V ■ {a{x, t)u{x, f)) -|- V ■ {d{x, t) ■ Vu{x, t)) -1- r{u, x, t), 
a; e D, f e [ 0 ,C]; a{x,t) = (0^(0;, t))j^i G d{x,t) = {dj{x,t)Yj^i G M.\ 

u{x, t) = gi{x, t), {x, t) G dQ X [0, C]; u{x, 0) = g2{x), x E Q, 

( 1 . 2 ) 

which is assumed to have some diffusion {dj{x,t) > do > 0, j = 
namely that it is not of pure hyperbolic type, and it is also assumed that some 
adequate spatial discretization based on Finite Differences or Finite Volume is 
applied to obtain the system fll.ip . Some stiffnes in the reaction part r{u, x, t) is 
also allowed. The treatment of Systems of PDEs do not involve more difficulty 
for our analysis but for simplicity of presentation we prefer to confine ourselves 
to the case of one PDE. 

We denote by Uhit) the solution of the PDE problem confined to the spatial 
grid (or well to the h-space related). It will be tacitly assumed that the PDE 
problem admits a smooth solution u{x,t) in the sense that continuous partial 
derivatives in all variables up to some order p exist and are continuous and 
uniform bounded on D x [0,f*] and that u{x,t) is continuous on D x [0,f*] 
(D = D1J(9D). It is also assumed that the spatial discretization errors 

ahit) := - fh(t, Uhit)), (1.3) 

satisfy in the norm considered, 

||c’‘h(Z:)|| (C > 0, r > 0), <t <t*, h —)■ 0. (1.4) 

In general C, C or C* will refer to some constants that maybe different at 
each occurrence but that all of them are independent of h —?■ 0 and from the 
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time-stepsize r —)■ 0. The vector norm used is arbitrary as long as it is defined 
for vectors of any dimension. For square matrices the norm used is the induced 
operator norm, ||y4|| = sup^^o ||Au||/||u|l. 

In spite of most of our results apply in general, we will provide specific results 
for weighted Euclidean norms of type 

It should be noted that in this case we have for any square matrix A that, 
||A|| = ll^lh, iV = l,2,3,.... 


We assume some natural splitting for fh (directional or other). 


fh{t,y) = (1.5) 

i=i 

which provides some natural splitting for the Jacobian matrix at the current 
point {tn,yn), 


Jt 


i=i 


Jh := 


^fh {iny 2/r? 

dy 


Jj,h ■- 


df j^hitny y-n) 
dy 


( 1 . 6 ) 


This goal of the paper is to analyze the convergence order of the Method of 
Lines (MoL) approach for time-dependent PDEs of Advection Reaction Dif¬ 
fusion PDEs, with the main focuss on the time integration of the large ODE 
systems resulting of the spatial PDE-semidiscretization, where some stiffness 
is assumed (parabolic dominant problems with stiff reaction terms) and the 
time integrators are based on very few iterations of splitting type (Approxi¬ 
mate Matrix Factorization and Newton-type schemes) applied to highly stable 
Implicit Runge-Kutta methods. It should be remarked that the underlying 
Implicit Runge-Kutta method is never solved up to convergence, hence the 
convergence study does not follows from the results collected in classical ref¬ 
erences about finite difference methods such as [miHllUlllTIITHII^ . The kind of 
approach to be considered here has interest since it is easily applicable to 
general systems of PDEs as we will see later on and it is reasonably cheap 
for non-linear problems in general (although we give convergence results for 
semilinear problems only) when some splitting of the function fh and its Ja¬ 
cobian is available and the split terms can be handled efficiently. In particular 
a method based on three AMF-iterations of the two-stage Radau IIA method 
[1] has shown to be competitive [7] when compared with some standard PDE- 
solvers such as VODPK [2]|3] in some interesting non-linear diffusion reaction 
problems widely considered in the literature. We also present two new meth¬ 
ods based on the 2-stage Radau IIA, by performing just one or two iterations 
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of splitting type, respectively. The method based on two iterations is one of 
the very few one-step methods of splitting type we have seen in the literature 
that has order three in PDE-sense for the time integration. 

The rest of the paper is organized as follows. In section 2 we introduce the 
AMFg-RK methods, and special attention is paid to some methods based on 
Radau IIA formulas. In section 3, the convergence for semilinear PDEs is 
studied in detail. The local and global errors are studied for the AMF^-RK 
splitting methods based on some general Runge-Kutta methods. Section 4 
is devoted to some applications of the convergence results to 2D and 3D- 
parabolic PDEs. 

Henceforth, for simplicity in the notations, we omit in many cases the h- 
dependence of some vectors such as fh, fj^h and of some matrices such as 
Jh and Jj^h U = l,...,d). It should be clear from the context which ones 
are h-dependent. Besides, we will refer to the identity matrix as I when its 
dimension is clear from the context. 


2 AMF-IRK methods 


For the integration of the ODEs fll.ll) . we consider as a hrst step an implicit s- 
stage Runge-Kutta method with a nonsingular coefficient matrix A = 
and a weight vector b = The method is given by the compact formula¬ 

tion (below <8) denotes the Kronecker product of matrices A®B = (aijB), A = 
(ciij), B = {pij)) 

Yn = e®yn A t{A (g) Im)F{Yn), 

Vn+l ^Vn “F (h (g IjyiyY^^ 

c = := Ae, e = (1,..., 1)^ G M®, := b^A~^, w = 1 — 

Yn = {Yn,j)j=i e M”"®, F{Yn) = (/(t„ + TCj, Kn,i))j=l ^ M”"®. 

( 2 . 1 ) 

It should be noted that we have replaced the usual formulation at the stepping 
point i/n+i = UnA T{b^ g Im)F{Yn) by the equivalent in fl2.1l) . which has some 
computational advantages for stiff problems when the algebraic system for the 
stages is not exactly solved. 

A typical Quasi-Newton iteration to solve the stage equations above is given 
by (below, J = df /dy {tn,yn) is the exact Jacobian at the step-point (tn,yn)), 

- .4 ® r J]A" = or', Y:^Y:-' + A\ v = (2.2) 
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where 


D:-^ ^ D{tn^ r, Vn, rr') ■■=e®yn- Y:-^ + t{A ® /™)F(Fr')- (2-3) 

A cheaper iteration of Newton-type when the matrix A has a multipoint spec¬ 
trum has been considered in piT^ (denoted as Single-Newton iteration) 

[Ims-n®Tj]A^ = D:-\ Y: = Y:-^ + A\ z/ = l,2,...,g (2.4) 

where 

T, = ^S,{I-L,)-^S-\ 7>0, 

Si, E are regular matrices and (2-5) 

Li, E M®’® are strictly lower triangular matrices. 

After some simple manipulations, by using standard properties of the Kro- 
necker product, this iteration can be rewritten in the equivalent form, 

[Is ® {Im - 1 tJ)]E'^ = ((A - L^)S-^ ® + {Lu ® Im)E'^, 

Yn = Yn~^ + {Si, (g) Iin)E’', u = l,2,...,q. 

To reduce the algebra cost, we use the Approximate Matrix Factorization [8] 
in short AMF, with J = and Jj = Jj^h given in fll.6l) . 

d 

:= Uilin - jrJj) = {Im - jtJ) + 0{t^), (2.7) 

i=i 

and replace in fl2.6p {I„i — l^J) by fl^i, which yields the AMF^-RK method 
based on the underlying Runge-Kutta method 


[h ® Ad]E^ = ((A - L,)S-^ 0 + {Lu ® Im)E'^, 

Yn = Yn~^ + {Si, 0 Im)E^i z/ = 1, 2,..., g 

Y^ = e®yn (Predictor) 
yn+i = ^yn + (B^ 0 Im)Y^ (Corrector). 


Our starting point for the convergence analysis in the next section takes into 
account that the AMF^-RK method can be rewritten in the equivalent form [3] 


[/ 0 / - T, 0 rP](F„" 
F;° = e 0 yn. 


Yn = D{tni r, yni Y^ l<v <q 
yn+i = wyn + {if 0 Im)Y^, 


(2.9) 
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where the matrix P plays a primary role 


P:= ( 7 r)-H/-n,) 

j<k j<k<l 

( 2 . 10 ) 


2.1 AMFq-RK methods based on the 2 stage Radau II A formula 


We are going to deserve special attention to AMF^-RK methods based on the 
2 stage Radau IIA formula [T]. This formula has coefficient Butcher tableau 
given by 



1/3 

5/12 

- 1/12 

1 

3/4 

1/4 


3/4 

1/4 


This is a collocation method (stage order is two) possessing good stability 
properties, such as L-stability (i.e. A-stability plus R{oo) = 0, with R{z) being 
the linear stability function of the method), and has order of convergence three 
(in ODE sense), not only on non-stiff problems but also in many kinds of stiff 
problems [1]. These properties for the underlying Runge-Kutta method are 
convenient, since the family of ODEs fll.ip involves stiffness in most of cases, 
due to the diffusion terms and possibly to the reaction part, and it is expected 
that the methods to be built on inherit part of the good properties of the 
original Runge-Kutta method. 


The next three AMF^-Rad methods have coefficient matrices {L^, and 
and eigenvalue 7 of the form 


T. = jSYl2-LY-^Sf\ S, 



. n = 

0 

0 

(01 j 


oy 


, 7 = Y^det(A) = l/\/ 6 . 


( 2 . 11 ) 

AMFi-Rad was derived in [5] by looking for good stability properties and order 
two (ODE sense). In particular the method is A( 7 r/ 2 )-stable for a 2-sphtting 
(see in Dehnition [T] below, the concept of stability for a d-splitting), A(0)- 
stable for any d-splitting and has stability wedges close to 6d = ti/{ 2{d — 1 )) 
for d = 3,4. The method is based on one iteration [q = 1) and was required 
to fulhl (A — Ti)c = 0 and it has coefficients given by 


Si = 


3 + 2^/6 

9 


h = ^(-12 + 5\/6). 


( 2 . 12 ) 
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AMF 2 -Rad was derived in [5] by looking for good stability properties and order 
three (ODE sense). The method is A( 7 r/ 2 )-stable for a 2-sphtting, A(0)-stable 
for any d-splitting and A( 7 r/ 6 )-stable for d = 3,4. The method is based on 
two iterations {q = 2 ) and their matrices Ti and T 2 were required to satisfy 
(A—Ti)c = 0 and {A—T 2 ) = 0^, = (0,1), respectively. Its coefficients 

are uniquely given by 


Si 


3 + 2x/6 


9 

5 


h = -(-12 + 5\/6) 


S2 = 


2\/6 


U = 


3^^Q 


(2.13) 


AMFs-Rad was derived in [laz] by looking for good stability properties and 
order three (ODE sense). The method is A( 7 r/ 2 )-stable for a 2-sphtting, A(0)- 
stable for any d-splitting and close to A( 6 *d)-stable for d = 3,4 with 6^ = 
7r/{2{d — 1)). The method is based on three iterations {q = 3) and their 
matrices T = Ti = T 2 = T 3 were required to satisfy e^T~^{A — T) = 0^. Its 
coefficients are uniquely given by 


5 -2x/6 , 

Si — S2 — S3 —---, li 


h = h 


3\/6 


(2.14) 


In [7], a variable-stepsize integrator based on the AMF 3 -Rad method was suc¬ 
cessfully tested on several interesting 2D and 3D advection diffusion reaction 
PDEs by exhibiting good performances in comparison with state-of-the-art 
codes like VODPK [2113] and RKC [T 6 |IT^ and its implicit-explicit counterpart, 
iRKc pms] . The other two methods, AMF^-Rad {q = 1,2), were introduced 
later [5] after carefully analyzing the PDE errors on semilinear problems and 
with the purpose of reducing the number of iterations w.r.t. AMF 3 -Rad. 


3 Convergence for semilinear problems 


For our convergence analysis we consider AMF^-RK methods applied to the 
ODE problems coming from the spatial discretizations of semilinear PDE 
problems of type fll. 2 p where the advection and diffusion vectors a{x,t) and 
d{x,t) are both constant and the reaction part has the form 

r{u, x,t) = K u + g{x,t), k being a constant, x G D C Mb (3.1) 
In this way, the ODE systems have the form 

y'hit) = fh{t, Vh) ■■= JhVhit) + gh{t), |/fe(0) = h 0+, 

Jh = Sj=l Jj,hi t G [0,t*]. 
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Here, the exact solution of the PDE conhned to the spatial grid Uh(t) = u{x, t), 
is assumed to satisfy fll.Sp and fll.dp . Thus, we focus on the global errors of 
the MoL approach, where the spatial discretization is carried out hrst by 
using hnite differences (or hnite volumes) and then the time discretization 
is performed by using AMFg-RK methods. It is important to remark that we 
will not pursue the details of the spatial semidiscretizations but rather it is 
assumed that the spatial semidiscretizations are stable and provides spatial 
discretization errors satisfying fll.dp . We shall provide uniform bounds for the 
global errors of the MoL approach {yh{t) henceforth denotes the numerical 
solution of the MoL approach) in the sense 

en,h-=Uh{tn)-yh{tn) = 0{TY^ + 0{h°‘T^^), h-)■ 0+, T 0+, (3.3) 

which is meant that there exist constants Ci, C 21 Pi, P 2 ; « (all of them inde¬ 
pendent on h and r) so that in the norm considered, 

||en,h|| < + C 2 h^TP\ h ^ 0+, r ^ 0+ holds. 

In our convergence analysis we need that all the matrices Jj^h pairwise conmute 
and that they can be brought to the following decomposition (it has some 
resemblance with the Jordan’s decomposition, but it is a little more general) 


Jj^h = Cond(0,j) := ||0,,|| ■ ||0;,^|| < C, h 0+, 1 < j < d, 

= BlockDia.g(A;y A'^ .... A<J>). A® = A® / + B®, 


Re A® < 0, 


dim(Ef) < N, \\E, 


f||oo<C", / = 1,2,...,79, (h^0+). 




are all of them strictly lower triangular matrices. 


(3,4) 


Another important approach for the convergence analysis of the MoL method 
(mainly concerned with the time integration) is based on the pseudo-spectra 
analysis of the matrix Jh [I3] and the related matrices Jj^h- That analysis is 
of more general scope but it is much more difficult to make and as we will see 
below, our analysis is enough for some interesting kind of semilinear problems 
and it is expected that the results extend to most of the non-linear problems 
of parabolic dominant type. 


Next, we consider a standard 3D semilinear-PDEs problem where the assump¬ 
tions in fl3.4p are fulhlled. 


3.1 An example 


Consider the semilinear PDE-problem fll.2l) with x G D = (0,1)^, with con¬ 
stant vectors, a{x,t) = d{x,t) = {dj)^=i, dj > 0 (j = 1,2,3) and 











r{x,u,t) as in fl3.ip . Consider the spatial semidiscretization by using second 
order central differences and spatial resolution h = 1/{N + 1). This yields a 
semilinear ODE systems of dimension m = of the form fl3.2p for d = 3. 
The matrices Jj^h are given by 


IN ^ -'TV 


Ji,h — In ® In J2,h — In In, J-i,h — T3 

Ti = Tridiag(«z, 61 , A) e / = 1 , 2 , 3, 

ai = h~‘^{di — 2~^h ai), / 3 ; = h~‘^{di + 2~^h ai), 61 = h~‘^{—2di + 

(3.5) 

and the vector ghit) includes the reaction part g{x,t) plus the boundary con¬ 
ditions. It is straightforward to see that the Ji^h pairwise commute. Moreover, 
by assuming Cell-Peclet numbers P p. 67, formula (3.42) ] 


h\ai\/di<2, I = 1,2,3, 

from [m section 2 ] it follows that their spectral decomposition has the form 


Ti = Tndiiig{ai,Si,^i) = ViAiVi, Vi = DiU, I = 1,2,3, 


U = 


Ai = I)iag{Xi^k)k=i > = 5i + 2^J^i cos 


kn 


kjTT \ 

N+U ^ 


T=i,iv 


D, = 


'TV+l 


iV+ 1’ 

is an orthogonal matrix and 

N 




(3.6) 


From here we conclude that all the matrices can be brought to the spectral 
decomposition in fl3.4p having negative eigenvalues and with matrix Qh = 
V 3 (D V 2 0 El. Observe that 


I| 0 hii 2 i| 07 'll 2 = nL llEibiif^-^ib = nL llAihiiA-'lh 

_ jlj- / 2 di -F ^^( 2 di + 

1=1 


\\ 2 di-h\ai\) l}\ 2 di-h\ai\) 

as h —)■ 0 . 


I 


3.2 Analysis of the Truncation Errors 


The AMFg-RK method applied on problem fll.ip can be expressed in the simple 
one-step format yn+i = 4>f{tn, yn, t), n >0. Thus, the time-space global errors 
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Uhitn) - Vn satisfy 


— 


(^n+l) 0/ (^n 5 (j'n) •) ^) ) ~l~ (0/ (j'n) •> 0/ (^n? Vu') ^) ) 

litw) 1 “t" [^^//*92/]n(^/i(^n) Vn)') 

where 

[50//5|/]„ = + (^ “ 'r)c?6', 

and the time-space local errors are defined by 

In — l(tni • ^/i(^n+l) ^h(tn)i (2''^) 

Then, we have for the time-space global errors the recurrence 

e^+i = [50//52/]n • Cn +/n, n = 0,1, 2,..., “ 1- (3-8) 

In order to get a better understanding of the latter recurrence, we next intro¬ 
duce the following matrix operators (P is dehned in fl2.10p l 

Qy = {I ® I - T^® tP)~^, = Qy{A ®tJ-T^® tP), z / > 1 , Qq = L 

(3.9) 


Lemma 1 The time-space global errors provided by the AMFg-RK method 
when applied to the problem A3.^) satisfy the recurrence 


en+l = Rq{Tj,TP) ■ Cn + ln, U = 0, 1, 2, . . . , P/t - 1, (3.10) 


where In stands for the time-space local error defined in \3. Tp and 


1 j 


RqirJ, tP) =wI + (B^ ® /) Q, + E(n (e 0 /), ( 3 . 11 ) 

V j=q i=q / 

with (vPl) given by ii3.9\) . Moreover, the function RqijJ^TP) fulfils 


Rq{Tj,TP) 


/ = (fi^ 0 J) + ^(n M,)Q,-1 -l[Mi 


j=q i=q 


i=q 


{c®tJ). (3.12) 


Remark 1 It must be observed that commutativity does not hold in general, 
thus Y\^j=q Mj = MqMq_i ■ ■ ■ Mi. Ou the other hand, Rq{-) can be seen as the 
linear stability function of the method. The identity fl3.12l) for the function 
Rq{-) — I will play a major role in a favourable propagation of the local errors 
in a similar way as indicated in Lemma 2.3 in [9l p.l62]. 


Proof of Lemma [T]. Our hrst step is to analyze the operator [dcff/dyjn for 
the semilinear problem fl3.2p . Taking into account that the method is dehned 
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by fl2.9p . then we are led to compute 


dyn+i 
dyn 


with yn+i = wyn + At 


this end, by taking derivatives with regard to yn in the iteration fl2.9p . it holds 
that 


{I ® I-T, 




= e®I + {-I ®I + A®tJ) 


dY 


u-l 


dyr. 


From here, after some simple manipulations it follows that. 


dY" dY"~^ dY^ 

^ = Q.(e<8)/) + M,^^, {u = l,2,...,q), -^ = e®I. (3.13) 

^Vn ^Vn ^Vn 

From an inductive argument, it is not difficult to see that 


dY^ 

dyn 


1 j 




j=q i=q 


Then, by denoting Rg{Tj,TP) := 


dyn+i 

dyn 


it follows 


Rg{Tj,TP) =zul + {if 


dY^ 

dyn 


(3.14) 


and we deduce both fl3.11l) and fl3.10l) from fl3.8p and fl3.14p . 


In order to prove fl3.12p . we first take into account that Rg{-)—I = (&" ®I)Z^, 
where Z" = dY"/dyn — e® I. Then, from the recurrence fl3.13p . it follows af¬ 
ter some simple calculations that Z'/ = MyZ'/~^+ Qy{c®TJ), (i/ = 1, 2,..., g), 

with = 0. From here, we deduce Z/^ = ( Q, + Edl - ri mA {c®tJ), 


J=q i=q 


i=q 


and this directly gives fl3.12p . 

Remark 2 For a given rational function of two complex variables 


□ 


Az.w) = 




mi 


m2 


-1 




— = E 




aijZ w 


v*j=o 


E 


(3.15) 


we define the associated mapping ({Z, W) for two arbitrary commuting matri¬ 
ces Z and W just by replacing z hy Z and why W whenever the denominator 
yields a regular matrix. Sometimes we are given the rational mapping (^(Z, W) 
first and then we define the rational complex function just by replacing the 
matrices Z and W by the complex variables z and w, respectively. The above 
definitions are straightforward extended to functions and mappings of more 
than two complex variables. 
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We will be mainly concerned with the case in which z = r J and w = tP, where 
J and P are dehned in fl3.2l) and fl2.10p . respectively. It should be noticed that 
for instance the (i, j)-element of the matrix M^, see fl3.9p . would be given by 
(observe that it is a matrix itself) 

MijijJ, tP) = {ej (g) I){Is ® Ira-Ty® tP)~^{A ®tJ -T^® ® 

where ej denotes the j-vector of the canonical basis in and the correspond¬ 
ing complex function is 

w) = eJ{Is - wT^)~^{zA - wT^)ej. 


Another important point is that despite of we are considering cases with a 
d-splitting for J as indicated in fl3.2l) . the replacement of every tJj by the 
complex variable Zj and the dehnition of 

Zk, tc := 7 "M 1 - 

fc=i v 

simplihes the study to the case of two complex variables 2 ; and w or well to 
the case of mappings acting on the two matrices tJ and tP. 

It is worth to mention that our rational mappings and related complex func¬ 
tions are all well dehned whenever Re 2 ;^ < 0 for k = 1, 2,..., d and d arbitrary, 
because the existence of the matrix inverse (J — wTy)~^ is guaranteed if and 

d 

only if (1 — jw)~^ ~ II exists. It is easily seen the existence of the 

k=l 

late expression by virtue of 7 > 0 and that all the eigenvalues of the matrices 
Jj, {j = 1,... ,d) have a non-positive real part. Moreover, for any u = 1,... ,q 
and any d > 1 , we next prove that 



11(1 -T'*) 


(3.16) 


fc=l 


sup \Qt,{z,w)\ <+ 00 , 

Re <0, k=l,...,d 


sup \M^{z,w)\ < -fcx). 

Re zic <0, k=l,...,d 


(3.17) 


2 ; and w dehned in fl3.16p . 

This see this, observe that (Tj,— 7 /) is a nilpotent matrix fulhlling {T^—'^iy = 
0 and that 


Q.,(z.w) = (I - wT.,) ^ = ((1 - wy)I - w(T., - yl)) ^ 

= (1 - »7)-' (/ -T^(Z- 7l)y‘ = (1 - »7)-‘ Ep; (t^)’ (r„ - yiy 

and 

AZ(z,w) = Q..(z, w)(zA - wZ) = (E-;J(^)l(r„ - 7/)^) 

- (E*i;(i^t(r. - 7/)i) r„. 
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Hence the boundedness of Qu{z,w) and Mu{z,w) follows from the bounded¬ 
ness of 


1 ^ 1 

= 1 n7i(i 

- 7^/ 

1 1—ii;7 1 

1 w 1 

1 

1 

1 1 

1 l—wy 1 

1—w'y 1 


\-ii 


< 7-^1 + !) = 27 


-1 


and from the next lemma. 


□ 


Lemma 2 For any d = 2,3,..and z and w defined in Hd.16]) . we have that 


sup 

Re zf^ <0 


1 — 7 W 


= 7 




d-i\ 1/2 


d^ 


:-2 


Proof. The third equality below follows from the Maximum Modulus princi¬ 
ple, which says that the Maximum Modulus is reached at the boundary of the 
open region for complex analytical functions, 


sup 

Re zi, <0 
k=l,...,d 


'yw 


= 7 ^ sup 

Re zz <0 

k=l,...,d 


'yz 


'yw 


= 7 ^ sup 

Re ui. <0 
k=i,...,d 


Ui U 2 Ud 


= 7 


-1 


{yi + y2 + ■■■ + yd)i 


nil v'l + (vkY- 


= 7 ^ max 

xi^ >0 

k=l,...,d 


uLii^-Uk) 

'' {xi + X2 + . . . + Xd)"^ 


uLia + ix^fi) 


The computation of the extrema by making zero the gradient of the real 
function of several variables (xi,... ,Xd) gives the maximum for xi = X 2 = 
... = Xd = {d — The proof follows after substituting above this value. 

□ 


Definition 1 A method of the form (12.9p is said to be A( 6 ')-stable for a d- 
splitting, if and only if 

\Rq{z, tc) I < 1, ^z, w given by fl3.16p whenever Zk G W( 6 '), k = 1,2,... ,d, 

where (we consider that the argument of a no-null complex number ranges in 

[-7r,7r)) 

W{9) := {m G C : M = 0 or |arg(—«)! < 6 }. (3.18) 


3.3 Analysis of the Local Errors 


Next, we study the time-space local errors In given by fl3.7p . We will see that the 
time-space local error is composed of two terms, l]^^ related to the predictor 
used in the AMFg-RK method and /W related to the quadrature associated 
with the underlying Runge-Kutta method. 
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Lemma 3 If the linear system has continuous derivatives ^ (t) up to order 
p + I in [0,t*] and the underlying RK method has stage order I >1 (i < p), 
i.e. 


Ad ^ = j , 


b^d-^=r\ J = l,2,...,i. 


Then, the local error R in (S/^) of the AMF^-RK method is given by 


L = + 

41 - (B’’ 0 /)(Q, + E]=,(ni.,Af.)Qj-i-n;.,Mi)£>„ + i„, (3.19) 

41 - (8^0/)(nL,M) Au,{t„), 


with 


Auhdr^ 


■= {Uh{tn + CiT) - Uh{tn))t=l = Ej=l ® I)Uh\d. 


rP+1 


pi 


[ {ci - ef+u^f’^^\tn + OT)de) 

JO J i=\ 


(3.20) 


and (we use, (a:)+ ■.= x ii x >t) and (a;)+ := 0 otherwise) 

Dn = ^ ^ {{d - jAd~^) ® u^hdn)) + f (^ip{9) 0 df’^^\tn + 9 t)) de+ 

j=i+i J- 

t{A ® I) {ah{tn + CiT))l^^ ; gj{0) = — {{ci- 0)% -p^aij{cj - 6')^"M 

A \ j=i /,=i 

^n= ^(1 - d)u^f\trf) + [ 4 >{ 6 ) + 67 ) 96 , 

(3.21) 

Proof. Let us define 

Dn ■= {Uh{tn+CiT))l^^-e®Uh{tn)-r{A®T){fh{tn+CiT,Uh{tn+CiT))l^^. (3.22) 

From fll.3p . it follows that 

Dn = {Uh{tn + CiT))l^^ - (m^( t^))-^^ - T(A 0/) (m), (t„ + Qt) - (T/^(t„ + Qt))-= 1 . 

(3.23) 

Now, by using the Taylor expansion with integral remainder (below C,{x) de¬ 
notes a generic function having r -F 1-continuous derivatives in an adequate 
interval) 

C(tn + X)=t^ 

1=0 


^r +1 /-I 

C^'\tn) + ^ / (1 - + 6 x)d6, (3.24) 

r! Jo 
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and applying it conveniently to Uhitn + Qr) and u'^^itn + Qr) in fl3.23p with 
r = p and r = p — 1 respectively, we deduce after some computations, the 
expression for £)„ in fl3.2ip . Observe that order stage i for the Runge-Kutta 
method implies that c^—jAc^~^ = 0, —1 = 0, j = 1,..., i. The expression 

for 6n is obtained in a similar way, but taking into account that this time we 
dehne, 

S 

dn := Uhitn + t) - WUhitn) “ + CjT). (3.25) 

j=l 

Let us now take Un ■= {uhitn + CiT))l^^ and := tJn — U^, where are the 
iterates obtained by the scheme (12.91) when the predictor = e ^ Uhitn) is 
taken on the exact solution of the PDE at tn, i.e. pn = Uhitn)- This gives as 
solution, see fl2.9p 

Vn+l = WUhitn) + (fi'^ ® I)Ul. (3.26) 

From (I3.25P and (I3.26p it follows 

In = Uhitn+l) - Vn+l = (h'^ O 1)^1 + K- (3.27) 

In order to compute we insert the expression for f/^ in (12.9p . It follows for 
the semi-linear problem (13.21) that 

(/®/-T^®rP)(A(;-A(;-1) =-T)(f„,r,M;,(f„),f/^-i) r to \ 

(z/ = l,2,...,g) 

= -(/ (g) / - A (g) r J)A(;-i + Dn, 

This implies that A(( = + QyDn, 1 < V < q, and from this recurrence 

a; = («,+E(ri - n M) £>„+(n Mi) a;, 

\ 3=q i=q i=q / *=<? 

with A° = Auhitn) in (13.201) . Now, from this expression and from (13.271) the 
formula fl3.19p follows. □ 

Theorem 1 Consider a family of matrices {Jk,h\k=i Ph, h —)■ O'*", as 
given in H3.S\) and h2.10i) . respectively. Assume that ifR^D holds and that 

d 

U Spect(Jfc,;,) c Wie), (h ^ 0+) (3.28) 

A:=l 

is fulfilled for some 9 G [0,7r/2]. Let Liz,w) he a complex rational function 
satisfying 


sup \Liz,w)\ <1, 2 ; and w given by (I3.16p . 

^fcGVV(0), k=l,2,...,d 

Then, we have that 

||L(r J, rP)"|| < (7*, 0 < nr < t*, (r, h —)■ 0+). 
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Proof. For simplicity of notations, we omit the sub-index h in the matrices. 
By virtue of fl3.2p . fl3.4p and fl3.15p it follows that 


II (L(r J, TP)r II = ||0 ■ (L(rA, rT))" ■ ©"i < CH (L(rA, rT))" |1, n > 1, 
where 

:= Yi=i rAk, Ak = Block-Diag(Ai^\ Ai^\ ..., A^^), 
tT ;=7-' {l-ntiil-^frAk)). 

By dehning rA^^^ := J2k=i and := 7 “^ (^I — 0^=1 (-^ — 7'7'Afc^)), for 
the norm considered it follows that 

II (L(rA, rT))" II = ^rnax^ II (^L(rA^^\ ||, n > 1. 

Consider any diagonal block A^^^ = + E {E = E^^^ for simplicity of 

notation. Observe that all the matrices E are strictly lower triangular and 
they have uniform bounded entries and uniform bounded dimensions, hence 
all of them are nilpotent with nilpotency index < N) and dehne 

d / d 

Zk = t\ 1 \ l<k<d, z = '^Zk, w = 7 "M 1 - n(l - 7 ^fc) 

k=\ \ ^=1 



it follows that, 

L(rA«, tTW) = L ij^{zkl + tE),^-\I - f[{I - i{zj + tE)) 

V/c^l k=l 

By dehning the function of d complex variables, 

/ d d 

-0(^1, . . . ,Wd) ■■= L 7“^(1 - n (1 “ 'y^k)) 

\k=l k=l 

we get that L{tA^’-'>= 'ip{zil + tE, ..., zj + tE). Then, by using the 
Taylor expansion for 'ip around r = 0 and taking into the nilpotency of the 
matrix E, we deduce that. 




tpizil tE, ...,ZdI + tE) = ^jJ{zl, 


5 Zd)I+ 

{Z1,Z2,. 


. d^dZd 


:Zd)- 


Now, since L{z,w) = L{zi,... ,Zd) and all its partial derivatives up to order 
N are uniformly bounded on the wedge yV{6), we can write that 


iP{ziI+tE, Zdl+rE) = tpi^zi, ..., ZdPl+rLp^,^, \\Lp^h\\ < C% {r,h^ 0+). 
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From here we get for 0 < rn < t* that 

II i;ipizil + tE, ... ,ZdI + tE))"' II = II (^L{z, w)I + || < (1+rC*)” < exp(rC*). 


□ 


3.4 Some mappings and definitions 


For a given mapping ^ C™’’™' where X and Y are two arbitrary square 

complex matrices of order m we dehne some associated mappings in the fol¬ 
lowing way, 


CW(X,F):=(C(X,F)-C(X,X))(F-X)-S whenever det(F - X) ^ 0, 
X) := lime^owhenever the limit exists. 

(3.29) 

In a recursive form, when det(y — X) 7 ^ 0 and C[*](X, X) exists, we continue 
by dehning 


CM(x,y):= (cW(x,y)-cM(V.>f))(y-^)-‘. 

CM(.Y,.Y):=lim„(,C''+"(A',A' + £/), ! = 1,2,..., i*. 

By assuming det(y — X) 7 ^ 0 and the existence of ^ = 1, 2,..., /*, 

it is straightforward to show by induction that 

C(x, X) = i: X){Y - X)' + y){y - Xf^'- (3.31) 

1=0 


We have considered for convenience that C^°](X, X) := ({X,Y). It should be 
noted that the commutativity of the matrices X and X is neither necessary 
in the dehnitions above nor in the formula fl3.3ip . 

To have a practical meaning of the mapping (CW(X, X) we show next that 
assuming Ci^yV) has I* continuous partial derivatives regarding the second 
variable, then it holds that 

CW(^, X), / = 1 , 2 ,..., r. (3.32) 

To see fl3.32p . we use the induction. For Z = 0 it is true for convenience. For 
Z = 1 it is true since 

CW(Y, Y) = limCl'kY, Y+e/) = limf-* (C(Y, Y + el) - C(Y, Y)) = |1(Y, Y). 
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Assume it is true up to /, we show it for / + 1 by using fl3.3ip in the second 
equality and the induction in the third equality below. The L’Hospital formula 
for limits (for the indetermination 0/0) is used /+! times in the fourth equality, 


A-) ^ A- + ./) ^ ,.m C(A.A + e/)-S| C^l(A.A)(./y 


= lim 
= lim 


C(A, A + £/) - A) 


d^^^C,{x,y) 


(/ + !)! 


(X,X + e/) = 




(Z + 1)! dy^^^ 


(A, A) 


These results can be trivially extended to vectors (and matrices), namely 
{Cij{X,Y)) G by applying them to each component (ij{X,Y) G 

C™’™. Sometimes we will make use of this kind of vectors as we will see in the 
next section. 


3.5 Bounds for the local errors 


The forthcoming convergence results for AMF^-RK methods are based in the 
Lemma II.2.3 P p. 162], which can be stated as follows 

Lemma 4 Assume that the global errors e„ = e„(r; h), of a one-step method 
satisfy the recursion ii3.10\) . where the local errors In can be split (uniformly 
on h and r) as 

In = {Rq{TJ,TP) — I) (j){tn)T^h°'TO{T'^h^), U = 0, 1, . . . , f*/r — 1, (3.33) 

where the function (j){t) and its first derivative regarding t are uniformly bounded, 
then the stability condition 


sup \\Rq{T.J,TPY\\<C, (3.34) 

l<n<t */t 
T->- 0 +, h^ 0 + 

implies that the global errors uniformly fulfil 

Cn = 0{T^h°‘)0{T’^h^), n = 1,... ,t*/ t, r —)■ O’*", h —)■ O'*". (3.35) 


General Assumptions on the semilinear problem. 

To bound the local errors and conseguently the global errors we henceforth 
assume that the exact PDE solution Uh(t) confined to the spatial grid and 
the semilinear problem liS.2\} fulfil / [i.51) -( [i'^[ ), i3-4\ l \3.2^) for some 6 G 

[0,7r/2], and that the following hypotheses (related the matrices J and P) hold 
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for some constants (not necessarily positive) ai, /3i and rj and some nonnega¬ 
tive integer I*, whenever h ^ and r —)■ O’*", 


Up - ^ r‘h» Oil), 1 ,., 

I (P — J C>(1) ^ .^ 

(P2) J^ul’(\t) = 0(1), fc = 1, 2,... + 1, for some 77 . 


(3.36) 


It should be noticed that always ao = 0, because the derivatives (up to some 
order) of the exact solution are uniformly bounded, i.e. u^h\t) = 0{1), t G 

[0,r], fc = 0,l,...,p+l. 

Theorem 2 Assume that the Runge-Kutta method has stage order £ and that 
sup \z/{Rq{z,w) — 1)\ < C, z and w given bv (13.161) . (3.37) 

2fcew(0), 

/c=l,2,...,d 

Then for the AMF^-RK method we have that, 

if" = 0(r//) + 0 (t'+‘), (t ^ 0, ft ^ 0), 

and 

/W = Th^OiX) + /+^(i?g(r J, tP) - I) (C>(1) + Th^^O{l)) , T^O, h^O. 


Proof. According to Lemma [3] the term /W of the local error is given by, 

l^^^=XrJ,rP)Dn + dn, (3.38) 

where 

/ 1 i 1 

XtJ, tP) := (6^®/) g,(r J, rP) + ^(H M(r J, rP))Q,_i(r J, rP) - H rP) 

V j=q i=q i=q 

(3.39) 

From Remark [2] we have that (e^ denotes the j-vector of the canonical basis) 
sup \f{z,w)ej\ <C, (j = l,...,s), 2 ;, w given by (13A^. 

Re 2t.<0 
k=l,...,d 

From Theorem [1] this implies that 

max ||,^(r J, rP)(ej ® 1)11 < C", r —)■ O'*", h —)■ O'*". 

Then, from 03.211) in Lemma [3] the hrst bound for /W follows. 
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For the second bound, we separate in fl3.38p the r^+^-term from the others, 
take into account fl3.39p and Lemma [3l we get 


4'’ = ^ {arJ, rP) ((c'+i -{i+ l)Ac^) ® /) + (1 - + ®, 

where ® = + O^tK^). 


(3.40) 

Next, we dehne the mapping (assume that J is regular only to simplify the 
proof) 


t.(T J, tP) := tF)-/)-‘ (e(Tj. tP) ((c'+‘ - (t + l)Ac‘) 0 /) + (1 - 6 ’'c'+‘)/) . 

(3.41) 

By using the assumption fl3.37l) . the bounds in Remark [2] and Lemma [2], it is 
not very difficult to see that 


sup \v{z,z)\ < +CXD. sup \zv^^\z,w)\ < + 00 , 2 ; and w given by fl3.16p . 
z£W(e) ^fc6W(e) 

k=l,2,...,d 

(3.42) 

Then, from fl3.40l) it follows that, 

4'’ = ^ tP) - I) v{tJ, rP)nr'H^n) + ® 

= ^ {RirJ, tP) - I) tJ) + vW(tJ, tP){tP - rJ)) + ® 

= ® + ^ {R{rJ^ tP) - I) v{tJ, r J)wr'^(t„) 

+ ^ (P(rJ,rP) - /)r;W(rJ,rP)(rJ)(J-i(P - 

= ® + ^ {R{rJ, rP) - I) Oil) + ^ (P(r J, rP) - I) Oirh^-) □ 

(3.43) 

For the analysis of the local error term in fl3.19l) . we dehne the mappings 


%iTJ,TX) := (13^0 J)n]=,M,(rJ,rX) G 
C,(r J, tX) := (P,(r J, rX) - J)"' V'q(r J, tX) e 
and their associated vector complex functions 

'ipqiz.w) := ‘ff'Y{]=qiI - wTj)~^izA - wTj) e 

Cq{z,w) := iRqiz,w) - l)~'^^|Jq{z,w) e C^’L 


(3.44) 


(3.45) 


These mappings will play a mayor role in the proof of the convergence results. 
It must be remarked that whereas ||'^q(2;, tc )||2 is uniformly bounded when z 
and w are given by 03.161) . the vector (qiz,w) = Oiz~^) as z —)■ 0 due to the 
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fact that (see fl3.12p t 


Rq{z,w) 


1 = Qq{z, w) + Mi{z, w))Qj-i{z, w) - n 


J=q i=q 


i=q 


(3.46) 

Hence '>^) is not bounded in general for z and w given by (I3.16p . However, 
Cq{z, z) is uniformly bounded as long as Rq(z, z) — 1 ^ 0 for z E W(6*)\{0}. 


From fl3.19p . by using fl3.3ip . we deduce that, 


/pi 


(RqirJ, tP) - I) E;=o TJ){h ® (r(P - J)y)^h{tn) 

{RqirJ.rP) - I)Cp\Tj,TP){h ® {t{P - J))'*+i)A,(t„). 


Next, we provide some convergence results for different kind of AMFg-RK 
methods, which depends on the Runge-Kutta method on which the AMF^-RK 
is based on. We start with Theorem [3] that meets applications for DIRK meth¬ 
ods (Diagonally Implicit Runge-Kutta) and SIRK methods (Single Implicit 
Runge-Kutta) and then with Theorems IU 0 and [ 6 ] which meet applications in 
the AMF^-Rad methods presented in section two. Of course, the assumptions 
(P1)-(P2) will be always assumed for some integers /* > 0 , £ > 1 , p > 1 . 


Theorem 3 IfT^ = A, u = 1,... ,q, with the Runge-Kutta coefficient matrix 
A having unique eigenvalue 7 > 0 (with multiplicity s), then the local errors 

(L = A lyy fulfil 


411 = 0{tK) + r^+HP(r J, rP) - /)(0(1) + O^rh^fi), 


42] _ / = 0,1,..., /, 


> (r O’*", h —)■ O’*"). 


I = max{ 0 , min{g — 2 , /*}}. 


If the method is A{6)-stable for a d-splitting and h3.2^) holds, then for any 
/ = 0 , 1, ..., /, the global errors fulfil (whenever r —)■ O’*" and h —)■ that, 

en,h = 0(/i4+r^min{l,max{r,r2/i/^i}}C>(l)+C>(r2'+2h^*+i); n = 1,2 ,... ,t*/r. 


Proof. The expression of /W was seen in Theorem [21 In order to show the 
expression for /pi, we start by deducing from fl3.44p and (13.121) that 


Cqiz,w) = {Rq{z,w) - 1) W{{I-wA) ^AY{z-wy, 
Rq{z, w) — 1 = (^{{I — wA)~^A f {z — wfi — /) {zA — I)cz. 


(3.48) 
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1 ^ 

From fl3.32p we have that C}^\z,z) = ■ 77 Tr~^(^, ^)- From here and from fl3.48p 

^ /! aw'- 

it follows that 

cj;\z,z) = o, / = 

From fl3.47p by taking I as npper index, for any / = 0,1,we have that 
= {R,{tJ,tP) - I)Cj;+^KrJ,rP)iIs 0 (r(P - 

= P{P,{tJ,tP) - I) {Cj;+^KTJ,TP){P(^TJ)) iP 0 J-\P - jy+^)Ah{P) 

= P{P{tJ, tP) - /)0(1)(/. 0 J-\P - jy+^){TC 0 u'M + 

= r2'+2h/3*+i(i?(rJ,rP) - 1)0(1). 

To see the bonnd for the global errors we apply Lemma HI The bonnds for 
the local errors P have been obtained above (see also Theorem [2] for /W). The 
bonndedness of the powers of PqijJ^TP) as indicated in fl3.34p follows from 
Theorem [1] by taking into acconnt the A( 6 *)-stability of the method for the d- 
splitting and that 03.281) holds. Now from Lemma H] the proof is accomplished. 

□ 

Theorem 4 For AMF^-PK methods with 7 > 0 and satisfying {A — Tpc = 0, 
we have that 

42 ] = r4P(r J, rP) - I) (C>(1) + h^^Oil)) , (r ^ 0+, h ^ 0+). 

Additionally if the method is A{6)-stable for a d-splitting and H3.28^) holds, 
then for r —)■ O'*" and h —)■ O'*", the global errors fulfil 

^n,h = C>(h'')+/min{l,max{r, r2h^i}}C>(l)+r2 (C>(1) + h^iC>(l)) ; n = 1,2,.. .,t*/T. 

Proof. The expression of /W was seen in Theorem [2l In order to show the 
expression for I'^f from 03.47p by setting /* = 0 we get that (observe that 
(fq{z,z)c = 0 because {A — Ti)c = 0. This expression is used in the third 
equality below) 

421 = {Pq{Tj,TP) - I)Cq{Tj,Tj)APtn) 

+ {Pq{Tj,TP) - I)Cjy\Tj,TP){P 0 (r(P - J)))APtn) 

= {PqpJ, tP) - I)Cq{Tj, tJ) (tC 0 J + r2C>(l)) 

+ (P,(r J, rP) - I) tP){I 0 r J)) {I, 0 {J-\P - J))) (rO(l)) 

= (P(r J, rP) - J) (r20(l)) + (P(r J, rP) - (r^h^^Oil)) . 
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This provides the bound for the local errors . The boundedness of the powers 
of RqirJ^TP) as indicated in fl3.34p follows from Theorem [T] by taking into 
account the A( 6 ')-stability of the method for the d-splitting and that fl3.28p 
holds. Now, from the bounds for the local error and from Lemma 0] the proof 
follows. □ 

Theorem 5 For AMFq-RK methods with 7 > 0 and satisfying 

sup \\z~'^Cq{z,z )\\2 <+ 00 , 

Jie 2 < 0, Z^O 

with 1 ] given in (P2) we have that 

4^' = (i?(r J, tP) - I) (0{t^+^) + 0{T^h^^)) , (r ^ 0+, h ^ 0+). 

Additionally if the method is A{6)-stable for a d-splitting and lid. 28\) holds, 
then for r —)■ O’*" and h —)■ O’*", the global errors fulfil 

^n,h = C>(h^)+min{l, max{r, n = 1 , 2 ,..., t*/r. 

Proof. The expression of /W was seen in Theorem [2l In order to show the 
expression for from fl3.47l) by setting /* = 0 we get that 

= {Rq{Tj,TP)-I)UTj,Tj)/Au{tn) 

+ {Rq{Tj,TP) - I)CI^\tJ,tP){R 0 (r(P - J)))Ah{tn) 

= {RqirJ, rP) - /)C,(r J, rJ) (rO(l)) 

+ {RqirJ, tP) - I) {C?{rJ, tP){I 0 r J)) {R 0 {J-\P - J))) {tO{1)) 

= {RqirJ,rP) - I) {C,q{rJ,rJ){I®{rJ)-R) (/0 (rJ)4 {rO{l)) 

+ {R{rJ,rP)-F)0{1) (r^h^^O{l)) 

= RqirJ, rP) - I) (0(1)) (r^+iJ 0 J’^O(l)) 

+ iRirJ,rP) - I) (r‘^h^^Oil)) 

= {RirJ, rP) - I) (0(r^+4 + Oir^h^R) . 

This provides the bound for the local errors The rest of the proof follows 
as in the previous theorems. □ 

Theorem 6 For AMFq-RK methods with 7 > 0 and 

(A-Ti)c = 0, sup^,,<o,^^o lk"\g(A 2:)||2 < + 00 , 
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with fj given in (P2) and assuming (PI) for I* = 1, we have that 

4^] = (i?(r J, tP) - I) (C>(r2+'?) + + 0{T^h^^)) , (r ^ 0+, ^ 0+). 

Additionally if the method is A{9)-stable for a d-splitting and 113.28\) holds, 
then the global errors fulfil 

^n,h = Oih"^) + min{l, max{r, T‘^h^^}}0{T^) + + 0{T^h°‘^) + 0{T^h^^), 

n = 1, 2,, t* /t, (r —)■ O’*', h —)■ O'*"). 

Proof. In order to show the expression for from fl3.47p by setting I* = 1 
we get that 

= {R,{rJ,rP)-I) (C,(rJ,rJ)+CW(rJ,rJ)(/®r(P-J))) 

+ (P,(rJ, tP) - I)C,f\TJ, tP){Is ®t\P- Jf)Ahitn) 

= {Rg{^J,rP) - /) (Cq{rJ,TJ) +Cfl(rJ,rJ)(/(8)r(P - J))) {{re® I)u'f,{tn) + 
+ (P,(r J, rP) - rP)(J, ® r2(P - J)2)(rO(l)) 

= (P,(rJ,rP) -/) (r2C,(rJ,rJ)(9(l) +rCfl(rJ,rJ)(J® r(P- J)O(l))) 

+ {R,{tJ,tP)-I) (cf(rJ,rP)(/,<8)rJ)) (/, <8) r J-^P - J)2)(rO(l)) 

= (P,(rJ,rP) - /) (r2 (C,(r J, r J)(r J)-*?) (rV’?0(l)) + 0(r3h“4) 

+ (P,(r J, rP) - /) (0(1) rV-4P - 
= {R{tJ, tP) - I) (0(r2+'?) + 0(r3h“i) + 0(r^h^=)) . 

This provides the bound for the local errors The rest of the proof follows 
as in the previous theorems. □ 


4 Application of the convergence resnlts for Dirichlet Bonndary 
Conditions in parabolic problems 

Let us next consider the 2D semi-linear diffusion-reaction model {e is a positive 
constant) 

Ut = ^(mxx + Uyy ) + 9{x, V, t), {x, y) G (0, 1 )^ t e [0, 1 ], > 0, ( 4 . 1 ) 

with prescribed Dirichlet boundary conditions and an initial condition. The 
PDE is discretized on uniform spatial meshes (a:,, yj) = {ih,jh), h = N~^, 1 < 
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hj < — 1, where iV — 1 is the number of interior grid-points for each spatial 

variable. We shall assume that the exact solution of the PDE fl4.ip is regular 
enough when {x, y, t) E [0,1]^ x [0, t*]. Let us denote Uh{t) := with 

a row-wise ordering, where Uij(t) := u{xi,yj,t) for 0 < i,j < N. Then, by 
using second-order central differences, we obtain for the exact solution of fl4.ip 
on the grid a semi-discrete system fll.2p with dimension m = {N — 1)^ 

= eJuh{t) + gh{t) + o-h{t) -f eh~‘^ur^{t), (4.2) 

where 

J '.= Ji + J 2 , Ji = In-1 <8 J 2 = Bm-1 ® In-1 -, ,, 

(4.3) 

Bn -1 = h-‘^TriDiag{l, -2,1) G h = 1/N. 


Moreover, ghit) = {g{xi,yj,t))^jj^, \\ah{t)\\ 2 ,h = 0{h^) (0 < t < T), whereas 
contains the values of the exact solution on the boundary, i.e.. 


(t) = it) 0 Cl -f it) 0 ejv-i -F Cl 0 wSf it) + cn-i 0 it ), (4.4) 


with uf^\t) 

and - 

in 




[u, 


n 

Above, {ei,..., eAr_i} denotes the canonical basis 


For the proof of the convergence results we need the lemma 0 and the lemma O 
given below. These lemmas can be derived from the material in [HI pp. 96-300] 
(see from Lemma 6.1 to Lemma 6.5). Lemmas [5] and [6] supply sharp values for 
the constants a;, and g appearing in the P-assumptions of section 3. These 
constants together with the convergence theorems provide specihc orders of 
convergence of the MoL approach for several AMF^-RK methods, in particular 
for the AMFg-Rad methods presented in section 2. 

The norm considered here for vectors, is the weighed Fuclidean norm 




J=l\\2,h 


\ 


1 

iW 


N-1 

E 


■^^3 \ 


= h\\{v.,)r-i\\2, 


and for matrices the corresponding operator norm. 


Lemma 5 Assume that exact solution u{x, y, t) of the 2D-PDE problem ^.1\ ) 
has as many continuous partial derivatives as needed in the analysis in (x, y, t) 6 
[0,1]^ X [0,f*]. Then for k = 1,2 ... and uj < \ we have that, 


want'd) 


2,h 

2,h 


0{1), and moreover 
0(1), whenever Mry(t) = 0. 
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Lemma 6 Assume that exact solution u{x, y, t) of the 2D-PDE problem lO 
has as many continuous partial derivatives as needed in the analysis in (x, y, t) G 
[0,1]^ X [0,t*]. Then, for I = 0,1,... we have that, 


iP-jMHt) 


2,h 


Oir^h^^), 




0{T^h^‘), 


where 


and 





(4.5) 

ai ^ 1 

[ — max{0, 3 -F 4(/ — 2)}, 

( — max{0, 3 -F 4(/ — 1)}, 

if ug {t) = 0, 
otherwise. 

(4.6) 

A = | 

[ — max{0,1 -F 4(Z — 2)}, 

1^ — max{0,1 -F 4(/ — 1)}, 

if ug {t) = 0, 

otherwise. 

(4.7) 




□ 


We next give a convergence theorem for 2D-parabolic PDEs when the MoL ap¬ 
proach with AMFg-Rad methods in section 2 are applied to the time discretiza¬ 
tion. The results still hold for 3D-parabolic problems (even for dD-parabolic 
problems and d > 3) and Time-Independent Dirichlet boundary conditions, 
but the proof requires some extra length to be included here. 

Theorem 7 The global errors (GE) in the weighted Euclidean norm of the 
MoL approach for the 2D-PDE when the spatial semi-discretization is 

carried out with second order central differences and the time integration is per¬ 
formed with AMPq-RK methods, are given in Table\^ There, q = min{l, 
and is meant for 0 {t^) where y < 2.25 is any constant. 


(r -)> 0+, h ^ 0+) 

GE (Time-Indep.) 

GE (Time-Dep.) 

AMFi-Rad 

0(/i2) + 0(r2) 

0{h^)+0{Q) 

AMF2-Rad 

Oih^) + OiT^) + T^O(g) 

Oih^)+0{Q) 

AMFa-Rad 

C)(/i2) +C)(r2-25‘) 

0{h^) + 0{Q) 


Table 1 

Global error estimates in the weighted Euclidean norm for Time-Dependent Dirichlet boundary conditions 
(in short Time-Dep.) and Time-Independent Dirichlet boundary conditions (in short Time-Indep.). 


Proof. In all cases we have that the stage order of the underlying Runge-Kutta 
Radau IIA method is i = 2 and the order of the spatial semi-discretization is 
r = 2. Moreover, all the three methods AMF^-Rad (g = 1, 2, 3) are A(7r/2)- 
stable for a 2-sphtting as it is shown in [5] for the cases q = 1 and q = 2 and 
in [7] for the case q = 3. Also, it should be noticed that fl3.28p holds. 

We start with the AM Fi-Rad method. We have for the case of Time-Independent 
Dirichlet Boundary conditions that the derivative regarding t vanishes on 
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boundary points {x^y) G F/i, i.e. = 0. From Lemma [6] we get that 

«! = 0 and /3i = 0. Then the bound for the global errors follows from Theo¬ 
rem m For the case of Time-Dependent Dirichlet Boundary conditions, from 
Lemma El we have that ai = —3 and /3i = —1. Then, the bonnd for the 
global errors follows from Theorem 01 The bonnd also applies to the AMF 2 - 
Rad method for Time-Dependent Dirichlet BCs, because this method fulhls 
the assumptions in Theorem 01 

For the case of the AMF 2 -Rad method and Time-Independent Dirichlet BCs 
we apply Theorem 0] for the case p = 1 and Theorem El with /* = 1 for the 
case g = Observe that from Lemma El we have that ai = 0 and /9i = 0 

and (32 = —1. Moreover the AMF 2 -Rad method fnlhls all the assnmptions in 
Theorem El by taking 1 ] = 1, see also Lemma El 

For the case of the AM Fa-Rad method and Time-Independent Dirichlet BCs 
we apply TheoremElwith any rj < 1.25, see LemmaEl Observe that in this case 
tti = 0, /9i = 0. Then from Theorem El the global errors are of size 0{h‘^) -F 
The proof that the order can be increased up to Oih?) + 
requires some extra technical details that we have omitted for simplicity. The 
case of Time-Dependent Dirichlet BCs follows from Theorem El too, but in 
this case (3i = —1. □ 


4-.1 Numerical Experiments 


We have performed some nnmerical experiments on two 2D-PDE and 3D- 
PDE problems of parabolic type in order to illustrate the convergence results 
presented in former sections for the AMF^-Rad methods. 

(1) Problem 1 is the 2D-PDE problem fl4.ip with diffusion parameter e = 0.1 
and Dirichlet Bonndary Conditions and an Initial Condition so that 

u{x^ y, t) = 10x(l — x)y{l — y)e^ + (3e'^^~'^~^, (4.8) 

is the exact solution. The case (3 = 0 provides Time-Independent Bonnd¬ 
ary conditions and no spatial error {(Thit) = 0, due to the polynomial 
natnre of the exact solution). The case (3 = 1 provides Time-Dependent 
boundary conditions and spatial discretizations errors of order two. 

(2) Problem 2 is the 3D-PDE problem (14. 9 h with diffusion parameter e = 0.1 


t) = e Au{lt, t) + g{lt, t), 
t G [0,1], ^ = (x, y, z) G (0, If G 


(4.9) 
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and Dirichlet Boundary Conditions and an Initial Condition so that 

u{x, y, t) = 64a;(l — x)y{l — y)z{l — z)e^ + (4.10) 

is the exact solution. Again, the case /9 = 0 provides Time-Independent 
Boundary conditions and no spatial error and the case /5 7 ^ 0 provides 
Time-Dependent boundary conditions and spatial discretizations errors 
of order two. 

On the end-point of the time interval t* = 1, in the weighted Euclidean norm 
we have computed as specihed in fld.lip . the global errors e 2 (h,r) (|/met(C) 
denotes the numerical solution at t* by the method considered), the number 
of signihcant hgures of the global errors 52 (h, r) and the estimated order of 
the global errors p(h, r) as powers of h when r = r/h is kept constant and 
both r and h tend to zero. 

e2{h,T):=\\uh{t*)-y^^t{t*)\\2,h^ 52(h, r) = - log^o e2(h, r) 

p(h, r) = {52{h/2, r/2) - 52{h, t))/ log^o 2. 


In the Tables min] and m we have considered for each h the time-stepsize r = qh 
for the corresponding AMF^-Rad method (g = 1, 2 , 3), so that all the methods 
make use of the same number of /-evaluations and similar CPU times in 
the computations. In those tables we have displayed the number of signihcant 
hgures in the global errors S 2 {h, r) and in brackets the estimated orders p{h, r) 
of each method. 

From Theorem [71 the global errors are expected to be of size (observe that 
r/h is kept constant) where: 

(1) for the AMFi-Rad method, y = 2 if Time-Independent BCs are consid¬ 
ered and /i = 1 if Time-Dependent BCs are imposed. This nicely hts with 
the results displayed in Table [2] (Time-Independent BCs) and in Table [3] 
(Time-Dependent BCs) for the 2D-PDE problem. Moreover, the conver¬ 
gence order is still y = 2 in the 3D-PDE problem for Time-Independent 
BCs as it can be seen in Table jH 

(2) For the AMF 2 -Rad method, /i = 3 if Time-Independent BCs are consid¬ 
ered and /i = 1 if Time-Dependent BCs are imposed. This hts well with 
the results displayed in Table [2] (Time-Independent BCs) and in Table [3] 
(Time-Dependent BCs) for the 2D-PDE problem. Moreover, the conver¬ 
gence order is also y = 3 in the 3D-PDE problem for Time-Independent 
BCs as it can be observed in Table jH 

(3) For the AMFs-Rad method, y = 2.25* if Time-Independent BCs are con¬ 
sidered and /i = 1 if Time-Dependent BCs are imposed. This can be 
observed in Table [2] (Time-Independent BCs) and in Table [3] (Time- 
Dependent BCs) for the 2D-PDE problem. Moreover, the convergence 




order also approaches to /r = 2.3 in the 3D-PDE problem for Time- 
Independent BCs as shown in Table 01 


h 

AMFi-Rad (p) 

AMF2-Rad (p) 

AMFs-Rad (p) 

r/h = 1 

r/h = 2 

r/h = 3 

1 

24 

82 = 3.74 (2.03) 

82 = 4.94 (2.82) 

82 = 4.90 (3.56) 

1 

48 

<52 = 4.35 (2.03) 

82 = 5.79 (2.89) 

82 = 5.67 (2.42) 

1 

96 

<52 = 4.96 (1.99) 

82 = 6.66 (2.92) 

82 = 6.40 (2.36) 

1 

192 

<52 = 5.56 (1.99) 

82 = 7.54 (2.93) 

82 = 7.11 (2.29) 

1 

384 

82 = 6.16 (2.03) 

82 = 8.42 (2.96) 

82 = 7.80 (2.29) 

1 

768 

82 = 6.77 (-) 

82 = 9.31 (-) 

52 = 8.49 (-) 


Table 2 

Significant correct digits (/ 2 ,h,-hiorhn) for th® 2D-PDE problem with Time-Independent Dirichlet BCs (/3 = 0). 
In brackets the estimated orders of convergence (by halving both the spatial resolution h and the time- 
stepizes r and taking ratio r = r jK). 


h 

AMFi-Rad (p) 

AMF2-Rad (p) 

AMF3-Rad (p) 

rjh = 1 

r/h = 2 

r/h = 3 

1 

24 

82 = 3.02 (1.00) 

82 = 2.79 (0.76) 

52 = 2.52 (0.66) 

1 

48 

82 = 3.32 (0.97) 

82 = 3.02 (0.83) 

52 = 2.72 (0.76) 

1 

96 

82 = 3.61 (1.00) 

82 = 3.27 (0.90) 

52 = 2.95 (0.86) 

1 

192 

82 = 3.91 (1.00) 

82 = 3.54 (0.93) 

52 = 3.21 (0.91) 

1 

384 

82 = 4.21 (1.03) 

82 = 3.82 (0.97) 

52 = 3.48 (0.97) 

1 

768 

82 = 4.52 (-) 

52 =4.11 (-) 

52 = 3.77 (—) 


Table 3 

Significant correct digits (^ 2 ,h"riorm) for the 2D-PDE problem with Time-Dependent Dirichlet BCs (/3 = 1). 
In brackets the estimated orders of convergence (by halving both the spatial resolution h and the time- 
stepizes r and taking ratio r = t jK). 


h 

AMFi-Rad (p) 

AMF2-Rad (p) 

AMFs-Rad (p) 

r/h = 1 

rjh = 2 

r/h = ‘8 

1 

24 

52 = 3.40 (2.03) 

52 = 4.31 (2.96) 

82 = 4.53 (2.69) 

1 

48 

82 = 4.01 (2.03) 

52 = 5.20 (2.96) 

52 = 5.34 (2.59) 

1 

96 

52 = 4.62 (-) 

82 = 6.09 (-) 

52 = 6.12 (-) 


Table 4 

Significant correct digits (/ 2 ,h,"riorm) for the 3D-PDE problem with Time-Independent Dirichlet BCs (/S = 0). 
In brackets the estimated orders of convergence (by halving both the spatial resolution h and the time- 
stepizes r and taking ratio r = r/h). 

As a conclusion we can say that the convergence results presented in Theo¬ 
rem [7] seem to be sharp for 2D-parabolic problems and that they still hold for 
dD-parabolic problems {d > 2) when Time-Independent boundary conditions 
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are considered. The proof of this fact requires some additional work and is 
not presented here. On the other hand, the convergence results are very poor 
when Time-Dependent Boundary conditions are considered. However, in such 
a situation we have developed a very simple technique (Boundary Correction 
Technique) to recover the convergence order as if Time-Independent Bound¬ 
ary conditions were considered. The explanation of the Boundary Correction 
Technique and the proof of the convergence orders requires some extra length 
and will be the objective of another paper. 

It is also important to remark that although we have considered in Theorem [TJ 
second-order central differences for the spatial discretization, the convergence 
results also hold for most of the usual spatial discretizations as long as they 
are stable and consistent with order r > 1. Numerical experiments carried by 
the authors seem to indicate that the convergence results also hold for many 
classes of non-linear problems. 
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